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ABSTRACT 



The emergence of tilted bipolar active regions (ARs) and the dispersal of their 
flux, mediated via processes such as diffusion, differential rotation and meridional 
circulation is believed to be responsible for the reversal of the Sun's polar field. 
This process (commonly known as the Babcock-Leighton mechanism) is usually 
modeled as a near-surface, spatially distributed a-effect in kinematic mean-field 
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dynamo models. However, this formulation leads to a relationship between po- 
lar field strength and meridional flow speed which is opposite to that suggested 
by physical insight and predicted by surface flux-transport simulations. With 
this in mind, we present an improved double-ring algorithm for modeling the 
Babcock-Leighton mechanism based on AR eruption, within the framework of 
an axisymmetric dynamo model. Using surface flux-transport simulations we 
first show that an axisymmetric formulation - which is usually invoked in kine- 
matic dynamo models - can reasonably approximate the surface flux dynamics. 
Finally, we demonstrate that our treatment of the Babcock-Leighton mechanism 
through double-ring eruption leads to an inverse relationship between polar field 
strength and meridional flow speed as expected, reconciling the discrepancy be- 
tween surface flux-transport simulations and kinematic dynamo models. 

Subject headings: Sun: dynamo, Sun: interior, Sun: activity 

1. Introduction 

Currently, some of the best tools for understanding the solar magnetic cycle are ax- 
isymmetric kinematic dynamo models and surface flux-transport simulations. On the one 
hand kinematic dynamo models (which are usually based on an axisymmetric formulation), 
attempt to model the magnetic cycle self-consistently by using a prescribed meridional flow, 
differential rotation, turbulent diffusivity and poloidal source (see Sec. [5]). They have been 
successful in reproducing several of the characteristics of the solar cycle (see for example: 
Choudhuri, Schiissler & Dikpati 1995; Durney 1997; Dikpati & Charbonneau 1999; Covas et 
al. 2000; Nandy & Choudhuri 2001; Rempel 2006; Guerrero & de Gouveia Dal Pino 2007; 
Jouve & Brun 2007; Munoz-Jaramillo, Nandy & Martens 2009, MNM09 from here on; for 
more information about kinematic dynamo models see review by Charbonneau 2005). On 
the other hand, surface flux-transport simulations study the evolution of the photospheric 
magnetic field by integrating the induction equation using a prescribed meridional flow, dif- 
ferential rotation and turbulent diffusivity. There are two main differences between surface 
flux-transport simulations and kinematic dynamo models: in the former the computational 
domain is restricted to the surface (without imposing axisymmetry) and they are not self- 
excited, but driven by the deposition of active region (AR) bipolar pairs. This type of models 
has proved a successful tool for understanding surface dynamics on long timescales (see, for 
example, Mackay, Priest & Lockwood 2002; Wang, Lean & Sheeley 2002; Schrijver, De Rosa 
& Title 2002) and the evolution of coronal and interplanetary magnetic field (see for ex- 
ample Lean, Wang & Sheeley 2002; Yeates, Mackay & van Ballegooijen 2008). However, a 
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discrepancy exists between kinematic dynamo models and surface flux-transport simulations 
regarding the relationship between meridional flow amplitude and the strength of the polar 
field (Schrijver & Liu 2008; Hathaway & Rightmire 2010; Jiang et al. 2010). On the one 
hand kinematic dynamo models find that a stronger meridional flow results in stronger polar 
field (Dikpati, de Toma & Gilman 2008), on the other hand surface flux-transport simula- 
tions find an inverse relationship (Wang, Sheeley & Lean 2002; Jiang et al. 2010). In this 
work we improve upon the idea proposed by Durney (1997) and further elucidated by Nandy 
& Choudhuri (2001) of using axisymmetric ring doublets to model individual ARs. We show 
that this captures the surface dynamics better than the a-effect formulation and resolves the 
discrepancy between dynamo models and surface flux-transport simulations regarding the 
relationship between meridional flow speed and polar field strength. 

2. Evolution of the Axisymmetric Component of the Magnetic Field on 
Timescales Comparable to the Solar Cycle 

As mentioned before, kinematic dynamo models are usually based on an axisymmetric 
formulation and our model is not an exception. Given that here we introduce an improved ax- 
isymmetric double-ring algorithm for modeling AR eruptions (see below), but AR emergence 
is strictly a non- axisymmetric process, it is important to study the amount of information lost 
by averaging over the longitudinal dimension. We do this by performing surface transport 
simulations driven by a synthetic set of AR cycles based on Kitt Peak data using the model 
of Yeates, Mackay & van Ballegooijen (2007). We perform a regular surface flux-transport 
simulation in which the bipolar ARs are distributed all across the surface of the Sun (Case 
1) and another in which the same set of ARs is deposited at the same Carrington longitude 
while leaving other properties (time, tilt, latitude of emergence and flux) intact (Case 2). 
The difference between both simulations is clear from the top row of Fig. [TJ where we show a 
snapshot of the surface magnetic field at the peak of the cycle for Case 1 (Fig. [l]-a) and Case 
2 (Fig. [l}b). Obviously these cases have entirely different magnetic configurations at the 
time of deposition. However, when the magnetic field is averaged in longitude and stacked 
in time to create a magnetic synoptic map (also know as butterfly diagram; Figs.[l}c&[l}d), 
a careful examination shows that the results are essentially the same within a margin of 1% 
(Figs. [l]-e & f). The reason the simulations have identical outcomes is that the differen- 
tial rotation and the meridional flow are both independent of longitude in the simulations. 
Note that non-axisymmetry is essential for the evolution of the corona and interplanetary 
magnetic field. This result simply indicates that an axisymmetric representation of surface 
dynamics is a reasonable approximation if we are only concerned with the general properties 
of the magnetic field at the surface over solar cycle timescales in the context of dynamo 
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models. 



3. Modeling Individual Active Regions as Axisymmetric Double- Rings 

The initial implementation of the double-ring algorithm by Durney (1997) and Nandy 
& Choudhuri (2001) consisted in searching the bottom of the convection zone (CZ) for 
places in which the toroidal field exceeds a buoyant threshold and placing two axisymmetric 
rings of constant radial flux directly above them. This implementation had two important 
deficiencies: strong sensitivity to changes in grid resolution and the introduction of sharp 
discontinuities in the (/) component of the vector potential. The first necessary step to address 
these problems is a careful mathematical definition of the vector potential associated with 
each ring doublet, which ensures a continuous first derivative in the computational domain. 
We do so by building a separable function: 

A ar (r,9) = K A($)F(r)G(e), (1) 

where K is a constant we introduce to ensure super-critical solutions and A(Q>) defines the 
strength of the ring doublet. F(r) is defined as 

f r<R Q -R ar 



sin 



-(r - (R Q - R ar )) 



where R e = 6.96 x 10 8 m corresponds to the radius of the Sun and R ar = O.85i? represents 
the penetration depth of the AR. This depth is motivated from results indicating that the 
disconnection of an AR flux-tube happens deep down in the CZ (Longcope & Choudhuri 
2002). G{9) ) on the other hand, is easier to define in integral form: 

G{6) = — f[B-{6') + B + {6')\ sm(e')de\ (3) 

sm " JO 

where B+ (£>-) defines the positive (negative) ring: 

o 9 < e ar t f - f 

B±(0) = { ±^y [l + cos(f(^-^±f))] ^Tf-f <e<0 ar Tl + f • (4) 



■sin(6>) 

o 9 > e ar t f + f 

Here 9 ar is the co-latitude of emergence, A is the diameter of each polarity of the doublet , 
for which we use a fixed value of 6° (heliocentric degrees) and \ = arcsin[sin(7) sin(A ar )] is 
the latitudinal distance between the centers, which in turn depends on the angular distance 
between polarity centers A ar = 6° and the AR tilt angle 7; x is calculated using the spherical 
law of sines (see Fig. [2]-a for a diagram illustrating these quantities). Fig. [5]-b shows the 
axisymmetric signature of one of such axisymmetric ARs. 
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4. Recreating the Poloidal Field 



Given that the accumulated effect of all ARs is what regenerates the poloidal field, we 
need to specify an algorithm for AR eruption and decay in the context of the solar cycle. On 
each solar day of our simulation we randomly chose one of the latitudes with fields higher 
than a buoyancy threshold of B c = 5 x 10 4 Gauss at the bottom of the CZ (r = O.71i? ), 
and calculate the amount of magnetic flux present within its associated toroidal ring. The 
probability distribution we use is not uniform, but is restricted to observed active latitudes. 
We do this by making the probability function drop steadily to zero between 30° (-30°) and 
40° (-40°) in the northern (southern) hemisphere: 



P{9) oc 1 + erf 



9 - 0.3057T 
0.055tt 



1 — erf 



9 - 0.694tt 
0.055tt 



(5) 



We then calculate the corresponding AR tilt, using the local field strength B , the calculated 
flux $ an d the latitude of emergence A. For this we use the expression found by Fan, Fisher 
& McClymont (1994): 

7 oc<X 5/4 sm(A), (6) 

reducing the magnetic field of the toroidal ring from which the AR originates. In order to 
do this, we first estimate how much magnetic energy is present on a partial toroidal ring 
(after removing a chunk with the same angular size as the emerging AR). Given that this 
energy is smaller than the one calculated with a full ring, we set the value of the toroidal field 
such that the energy of a full toroidal ring filled with the new magnetic field strength is the 
same as the one calculated with the old magnitude for a partial ring. Finally, we deposit a 
double-ring (as defined in Section [3]) with these calculated properties, at the chosen latitude. 



5. The Kinematic Mean-Field Dynamo Model 



We perform dynamo simulations to explore how the double-ring formulation compares 
to the near surface a-effect formulation. In particular we focus on the relationship between 
meridional flow speed and polar field strength. Our model is based one the axisymmetric 
dynamo equations: 



OB 
~di~ 



v p - V 



OA 
~dt 

B 



+ - s [v p • V(sA)] = r] ^V 2 - -^J A + a /(r, 6)F{B tc )B tc 



+(V-v p )B = v V 



B+s([V x (Aet)] -Vtt)- 



1 d(sB) dr] 
s dr dr 



(7) 



, (8) 
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where A is the ^-component of the vector potential (from which B r and Bq can be obtained), 
B is the toroidal field {B ( f ) ) ) v p is the meridional flow, Q the differential rotation, rj the 
turbulent magnetic diffusivity and s = rsin(#). The second term on the right-hand side 
of Equation [7] corresponds to the poloidal source in the mean-field formulation. In this 
formulation a is a constant that sets the strength of the source term and is usually used to 
ensure super-critical solutions; a{r ) 9) captures the spatial properties of the BL mechanism: 
confinement to the surface, observed active latitudes and latitudinal dependence of tilt, while 
F(B tc ) adds nonlinearity to the dynamo by quenching the source term for values of the 
toroidal field at the bottom of the CZ B tc that are too strong or too weak. More information 
about this source can be found in MNM09. Note that for simulations using the double-ring 
algorithm this term is not present in the equations (a = 0). 

In order to integrate these equations, we need to prescribe four ingredients: meridional 
flow, differential rotation, the poloidal field regeneration mechanism, and turbulent magnetic 
diffusivity. For the differential rotation, we use the analytical form of Charbonneau et al. 
(1999), with a tachocline centered at O.7i? whose thickness is O.O5i? and we use the 
meridional flow profile defined in MNM09. This meridional flow better captures the features 
present in helioseismic data, specially the latitudinal dependence. We use an amplitude of 
20 m/s for the results shown in Figure [3] and a variable amplitude for the results shown in 
Figure [5] (see below). We use a double stepped diffusivity profile as described in MNM09. It 
starts with a diffusivity value rj^d = 10 8 cm 2 /s at the bottom of the CZ, jumps to a value 
of rj cz = 10 11 cm 2 /s in the CZ, and then to a value of rj sg = 10 12 cm 2 /s in the near-surface 
layers. The first step is centered at r cz = O.71i? and has a half-width of d cz = O.O15i? and 
the second step is centered at r sg = O.95i? and has a half- width of d sg = O.O25i? . 

For the poloidal field regeneration mechanism we use the improved ring-doublet algo- 
rithm described above, using a value of K = 400, in order to insure super-criticality (for a 
meridional flow of 30 m/s). For those simulations which use an a-effect formulation, we use 
the non-local poloidal source described above (more information in MNM09) using a value 
of a = 0.25, in order to insure super-criticality (for a meridional flow of 30 m/s). 

6. Addressing the Discrepancy Between Kinematic Dynamo Models and 

Surface Flux- Transport Simulations 

In order to have a net accumulation of unipolar field at the poles, it is necessary to 
have an equal amount of flux cancellation across the equator. Since the meridional flow is 
poleward in the top part of the convection zone, it essentially acts as a barrier against flux 
cancellation by sweeping both positive and negative AR polarities towards the poles resulting 
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in weak polar fields. This leads to an inverse correlation between flow speed and polar field 
strength which is accurately captured in surface flux transport simulations. Contrarily, 
dynamo simulations in typically used parameter regimes obtain an opposite relationship not 
consistent with the above physics. This is because if there is already a strong separation 
of flux, a fast meridional flow will lead to an enhancement of the polar field due to flux 
concentration. This unrealistically strong separation is typical of kinematic dynamo models 
that use a non-local a-effect BL source (see Fig. [3]-c). The reason is that by increasing the 
vector potential A proportionally to the toroidal field B at the bottom of the CZ (Eq. one 
creates strong gradients in the vector potential above the edges of the toroidal field belt; this 
ends up immediately producing poloidal field which is as large in length scale as the toroidal 
field itself, circumventing the whole process of flux transport by circulation and diffusion. 
Figure [3] illustrates this fundamental difference: The top row shows the evolution of the 
surface magnetic field for a dynamo model using the double-ring algorithm (Fig. [3]-a) versus 
one using the a-effect formulation (Fig. [3[b). The different way in which each formulation 
handles the surface dynamics is evident. The double-ring simulation clearly shows a mixture 
of polarities and small-scale features which migrate to the poles (very much like the observed 
evolution of the surface magnetic field). On the other hand, the mean field formulation only 
shows two large scale polarities whose centroids drift apart as the cycle progresses. The 
bottom row depicts a snapshot of the poloidal field for the double-ring algorithm (Fig. [3^ 
c) and the a-effect formulation (Fig. |3[d) - both snapshots taken at solar max. Here the 
presence of small-scale features and a mixture of polarities is evident for the double ring, 
whereas the a-effect formulation only shows a large-scale magnetic field with two polarities. 
It is clear that although the large scale internal field is similar for both, the double-ring 
algorithm does a much better job of capturing the surface dynamics. 



6.1. Polar Field Strength vs. Meridional Flow Speed 

In order to study the relationship between meridional flow and polar field strength, we 
perform simulations in which we randomly change the meridional flow amplitude from one 
sunspot cycle to another (between 15 — 30 m/s). This is illustrated in Fig. [i] where a series 
of sunspot cycles is plotted along with their associated meridional flow. We then evaluate 
the correlation between the amplitude of the meridional flow of a given cycle and the polar 
field strength B r at the end of it. Since we want to evaluate the relative performance of the 
double-ring algorithm as opposed to the non-local BL source, we perform the same simulation 
for both types of sources. Aside from the varying meridional flow amplitude and the poloidal 
source, the rest of the ingredients are the same. It is important to note that partly due to 
difficulties in tracking the exact occurrence of solar minimum, the two hemispheres eventually 
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drift out of phase in long simulations - sometimes this phase difference leads to quadrupolar 
solutions which often go back to the observed dipolar solution. This parity issue only appears 
when the meridional flow is changed at solar minimum: if there are no variations, or if the 
variation takes place at solar maximum, the cycle is always locked in phase with dipolar 
parity. Nevertheless, to compare our simulations with surface flux-transport models, we 
change the flow speed only at solar minimum. To be consistent, we accumulate statistics 
only from cycles in which the two hemispheres are in dipolar phase. The statistics performed 
for both types of source contain about 200 sunspot cycles. 

The values of polar field we find using the kinematic dynamo simulations are of the order 
of 10 kG which is a common feature of dynamo models, which are successful in simulating 
the strong toroidal field necessary to produce sunspots and sustain the solar cycle (Dikpati 
& Charbonneau 1999; Chatterjee, Nandy & Choudhuri 2004; Jiang & Wang 2007; Jouve at 
al. 2008). Recent high resolution observations of the polar region have now confirmed the 
existence of such strong kilo-Gauss unipolar flux tubes (Tsuneta et al. 2008). Fig.[5]shows the 
results of both simulations. We find a weak positive correlation between meridional flow and 
polar field strength for the simulations using the non-local a-effect formulation (Fig. [5]-top), 
which is in general agreement with the results of Dikpati, de Toma & Gilman (2008). On the 
other hand, the simulations using the double-ring formulation distinctively show a negative 
correlation (Fig. j5[bottom), in agreement with surface flux-transport simulations (Wang, 
Sheeley & Lean 2002). This clearly establishes that the discrepancy between the models is 
resolved by introducing the double-ring algorithm and that the double-ring formalism does 
a better job at capturing the observed surface magnetic field dynamics than the non-local 
a-effect formalism. 



7. Concluding Remarks 

In the first half of this work, we perform surface flux-transport simulations to test the 
validity of the axisymmetric formulation of the kinematic dynamo problem. Our results sug- 
gest that this axisymmetric formulation captures well the surface flux dynamics over spatial 
and temporal scales that are relevant for the solar cycle. Building upon this we introduce an 
improved version of the double-ring algorithm to model the Babcock-Leighton mechanism for 
poloidal field regeneration in axisymmetric, kinematic dynamo models. We show that this 
new double-ring formulation generates surface field evolution and polar field reversal which 
is in close agreement with observations. Additionally, we find that this improved treatment 
of the Babcock-Leighton process generates an inverse relationship between meridional flow 
speed and polar field strength - which is suggested by simple physical arguments and also 
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predicted by surface-flux transport simulations. This resolves the discrepancy between kine- 
matic dynamo models and surface flux-transport simulations regarding the dynamics of the 
surface magnetic field. Since the latter drives the evolution of the corona and the heliosphere, 
our work opens up the possibility of coupling dynamo models of the solar cycle with coronal 
and heliospheric field evolution models. 
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Fig. 1. — Long term evolution of the photospheric magnetic field in surface flux-transport 
simulations: (a) Snapshot of the magnetic field at the peak of the cycle for Case 1 (active 
regions deposited at all longitudes), (b) Snapshot of the magnetic field at the peak of the 
cycle for Case 2 (active regions deposited at a single longitude), (c) Butterfly diagram for 
Case 1. (d) Butterfly diagram for Case 2. (e) Difference between the butterfly diagrams of 
Case 1 and Case 2. (f) Longitudinal average of the snapshots shown in the top row, taken 
at the peak of the cycle. 
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Fig. 2. — (a) Diagram illustrating the quantities which define the latitudinal dependence 
of a double-ring bipolar pair, (b) Poloidal field lines of one of our double-rings including 
a potential field extrapolation for the region outside the Sun. The dashed line marks the 
location of the penetration depth R ar . 



Double-ring Algorithm 



a-effect Formulation 




Fig. 3. — Comparison between surface dynamics as captured by the double-ring algorithm 
(left column) and the a-effect formulation (right column) . The top row shows the evolution 
of the surface magnetic field in the form of synoptic maps - the colormap is saturated to 
enhance the visibility of the field at mid to low latitudes. The bottom row shows a snapshot 
of the poloidal components of the magnetic field taken at solar max. The solid contours 
corresponds to clockwise field-lines, the dashed contours correspond to counter-clockwise 
field-lines. The thick dashed lines mark the location of the tachocline. 
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Fig. 4. — Diagram for the evolution of the meridional flow amplitude with respect to the 
sunspot cycle: each solar cycle TV has a unique meridional flow strength v n which is randomly 
chosen between 15— 30 m/s. Additionally, the polar field strength B r of cycle TV it is measured 
at the end of it. 
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Fig. 5. — Relationship between randomly varying meridional flow speed and polar field 
strength for simulations using the mean-field formulation (top row) versus simulations using 
the double-ring algorithm (bottom-row) . The polar field strength (in Gauss) is represented by 
the maximum amplitude of the polar radial field (Br) attained during solar minimum. The 
relationship between the above parameters is determined by the Spearman's rank correlation 
coefficient. Top-row: (correlation coefficient, r 0.325, confidence, p 99.99%). Bottom-row: 
(r -0.625, p 99.99%) 



